Constellation Plots: Neocortex, 2nd Trimester

All Samples, All cell types

LOAD DATA.

______________________________________________________________________________

• Neocortex_v3.RData Not necessary if just plotting on already calculated PCA, UMAP, clusters.

# Full dataset:
load("../data/cfc4b2_neocortex_v3.RData")

# 40k random subset (toy dataset)
# ncx.40k <- read_rds("constellation_plots/ncx_v3_40k.rds")

# 100k random subset of excitatory neuron lineage only.
ncx.exn.100k <- read_rds("../data/exn_lineage/toy/ncx.v3.exn.sub_100k.rds")

• Cluster / marker tables.

ncx.clusters <- read_delim("../tbls/83d19_Neocortex_allindividuals_combinedclusters_v1.txt",
                              "\t", escape_double = FALSE, trim_ws = TRUE)

ncx.markers <- read_delim("../tbls/8f0fc_Neocortex_subset1_clustermarkers_combo2.txt", 
                  "\t", escape_double = FALSE, trim_ws = TRUE)

********************************************************************

1. START scrattch.hicat functions to build constellation plots.

********************************************************************

1. Define cells to build the plot from.

[All cells: 404.2K]

a. All clusters / cell types:

Keep only cells with combo2 $combined.cluster.2 annotation.

cells.cl.df <-  ncx.clusters %>% filter(str_detect(combined.cluster.2, "combo2"))

s.obj <- Neocortex %>% subsetSeurat(cells.keep = cells.cl.df$cell.name)
# 348K cells

# Sanity check
cells.cl.df$cell.name %>% identical(s.obj@meta.data %>% rownames)
[TRUE]

b. Excitatory neuron lineage only:

Keep only cells with combo2 $combined.cluster.2 annotation AND whose $combined.cluster.2 annotation belongs to excitatory neuronal lineage classes.

cells.cl.df <- ncx.clusters.exn <- 
                  ncx.clusters %>% 
                      filter(str_detect(combined.cluster.2, "combo2") ) %>% 
                        filter(str_detect(combined.cluster.2, "Neuron|CR|Dividing|RG|IPC|OPC") )
# 271.4K cells in excitatory lineage.

b.1. Cells in exn linage 100k cell toy object:

# Define which object the constellation plot is based on:
s.obj <- ncx.exn.100k

cells.cl.df <- ncx.clusters %>% 
                  filter(cell.name %in% rownames(s.obj@meta.data) )

2. Build dataframes for constellation plots.

1. cl

cl <- cells.cl.df$combined.cluster.2 %>% as.factor %>% 
          set_names(cells.cl.df$cell.name)
# 128 clusters in all _combo2_ cells
# 77 clusters in all _combo2_ & ExN lineage cells.

2. rd.dat

rd.dat <- list(umap = s.obj@reductions$umap@cell.embeddings,
               pca = s.obj@reductions$pca@cell.embeddings)

# Subset UMAP and PCA cell embeddings: keep only cells in ncx.clusters.
# (Exclude cells in cluster 0 and keep only cells with "combo2" in combined cluster 2 name.)
# rd.dat %<>% lapply(function(x) x[names(cl), ])

3. cl.df

cl.df <- get_cl_df(cl)

cl.df$clade <- str_split_fixed(cl.df$cluster_label, "_", 2)[ ,1] %>% tolower 
# Add clade_id, clade_color to cl.df
clade.cols <- data.frame(# clade = unique(cl.df$clade),
                        cluster_color = c("cr"= "darkgrey", 
                                  "dividing" = "darkkhaki", 
                                  "neuron" = "deepskyblue", 
                                  "inteneuron" = "deeppink2", 
                                  "ipc" = "brown4", 
                                  "microglia" = "darkorchid1",
                                  "opc" = "cadetblue", 
                                  "other" = "darkslateblue", 
                                  "rg" = "darkorange", 
                                  "vascular" = "blanchedalmond")
            ) %>% rownames_to_column("clade")

cl.df %<>% left_join(clade.cols)
rm(clade.cols)

# cells.cl.df: Add cluster_id column from cl.df; remove unused columns. 
cells.cl.df <- left_join(cells.cl.df %>% select(cell.name, cell.type, combined.cluster.2),
                           cl.df %>% select(cluster_label, cluster_id), 
                           by = c("combined.cluster.2" = "cluster_label")
                ) %>% mutate(cluster_id = as.factor(cluster_id))

4. rd.cl.center Find cluster centers from UMAP coordinates

rd.cl.center <- get_RD_cl_center(rd.dat$umap, cl) 

rd.cl.center %<>% 
  as.data.frame %>% 
  set_names(c("x", "y")) %>%
  add_column(cl = 1:nrow(rd.cl.center), .before = "x") %>%
  # add_column preserves rownames.
  # but moving rownames to column cluster_label anyway bc of left_join below.
  rownames_to_column("cluster_label")

5. cl.center.df

Join cl.df and rd.cl.center into cl.center.df for input into get_KNN_graph.

cl.center.df <- left_join(rd.cl.center, cl.df,
                           by = c("cluster_label")) 

6. knn.cl

# Fixes needed for proper output of knn.cl.df
# cl.center.df %<>% rename(cluster_size = "size")
# levels(cl) <- cl.df$cluster_id
cl.numeric <- cl
levels(cl.numeric) <- cl.df$cluster_id

knn.result <- RANN::nn2(data = rd.dat$pca[, 1:10], k = 15)

knn.cl <- get_knn_graph(rd.dat = rd.dat$pca[ , 1:10], 
                        cl.df =  cl.df, cl = cl.numeric, 
                        k = 15, 
                        knn.outlier.th = 2, outlier.frac.t = 0.5)

rm(rd.dat, ncx.clusters)

3. Make constellation plot

# Keep only cells where $frac [fraction of cells in cluster with nearest neighbors in different cluster] >= 0.05.
# Defined in `get_knn_graph`: 
# knn.cl.df$frac = knn.cl.df$Freq / knn.cl.df$cl.from.total
# 10% : 213 edges
knn.cl.df <- knn.cl$knn.cl.df  %>% filter(frac >= 0.1)

# Plot only edges between ExN lineage clusters.
# knn.cl.df %<>% filter_at(vars(cl.from.label, cl.to.label), 
#                        all_vars(str_detect(., "RG|IPC|Neuron|OPC|Dividing"))
#              )

cl.center.df$cluster_label %<>% str_remove("_combo2")
  
cl.plot <- plot_constellation(knn.cl.df, cl.center.df, out.dir = "../out",
                              node.label = "cluster_label", exxageration = 1, curved = TRUE, 
                              plot.parts = FALSE, plot.hull = NULL, 
                              plot.height = 40, plot.width = 40,
                              node.dodge = TRUE, label.size = 2, max_size = 20)

——————————————————————————

2. Find DE genes between connected clusters.

——————————————————————————

# Dataframe w/ the proportion of k(15) nearest neighbors in each cluster for every cell.
nn.cl.df <- knn.cl$pred.result$pred.prob %>% as.data.frame

cl.df$cluster_id %<>% as.factor()
# Possibly move to top, before making all DFs.
cells.cl.df %<>% rename(cluster_label = "combined.cluster.2") %>% 
                  mutate(cluster_label = str_remove(cluster_label, "_combo2"))

cl.df %<>% mutate(cluster_label = str_remove(cluster_label, "_combo2"))


# Add column with cells' own cluster assignment from `cells.cl.df`.
nn.cl.df %<>% left_join(cells.cl.df,
                        by = c("query" = "cell.name")
                        )
# Add cluster_label corresponding to nn.cl.

nn.cl.df %<>% left_join(cl.df %>% select(cluster_label, cluster_id),
                        by = c("nn.cl" = "cluster_id"))

nn.cl.df %<>% select(query, cluster_id_self = "cluster_id", 
                    cluster_label_self = "cluster_label.x",
                    cluster_id_nn = "nn.cl",
                    cluster_label_nn = "cluster_label.y",
                    freq = "Freq")

# nn.cl.df %<>% select(query, combined.cluster.2, cluster_id, nn.cl, )

x <- filter(knn.cl$knn.cl.df, frac >= 0.1 & cl.from != cl.to) %>% arrange(cl.from)

knn.cl$knn.cl.cl.counts %>% head

cl.df

x <- filter(nn.cl.df, cluster_label_self == "Neuron_8" & cluster_label_nn == "Neuron_3")
# 12,201 cells in Neuron_8

x %>% filter(freq > 0) %>% 
  ggplot() + geom_density(aes(freq))

Neuron_3 cells with nearest neighbors in Neuron_8

x <- cell.cl.counts %>% filter(cluster_label == "Neuron_3") %>% rename(neuron_8 = "44")

median.nnCounts <- median(x$neuron_8[x$neuron_8 > 0], na.rm = TRUE)

  ggplot(x) + 
    ggtitle("neuron_3 cells \n n/15 nearest neighbors in neuron_8") +
    
    geom_bar(aes(x = neuron_8)) +
    geom_vline(xintercept = median.nnCounts, colour = "red") +
    annotate("text", x = median.nnCounts + 1, y = 70,
             label = paste0("median = ", median.nnCounts)
  ) + 
    xlab("# of neuron_8 neighbors") +
    ylab("") +
    theme_minimal() 

NA

Compare cells above and below the median count of neuron_8 neighbors.

cells <- list()
cells$above.median <- x %>% filter(neuron_8 > median) %>% pull(cell.name)
cells$below.median <- x %>% filter(neuron_8 < median)  %>% pull(cell.name)

s.obj <- SetIdent(Neocortex, cells = cells$above.median, value = 'nn_above_med') %>% 
          SetIdent(cells = cells$below.median, value = 'nn_below_med')
# levels(s.obj)
markers <- list()
markers$nn_above_med <- FindMarkers(s.obj, slot = "scale.data", 
                                  # cells.1 = cells$above.median, cells.2 = cells$below.median,
                                  ident.1 = "nn_above_med", ident.2 = "nn_below_med", 
                                  logfc.threshold = 0)

write_tsv(markers$nn_above_med, "../out/DEgenes_neuron3_vs_neuron8.tsv")

Calculate enrichment ratio, filter genes w. adj p-value < 0.05, sort table.

markers.tmp <- markers$nn_above_med %>% rownames_to_column("gene") %>%
  mutate(enrich.ratio = pct.1 / pct.2,
         gene.score = avg_diff * enrich.ratio,
         across(.cols = where(is.numeric), .fns = round, digits = 5)
  ) %>%
  filter(p_val_adj <= 0.05) %>% 
  select(gene, pct.1, pct.2, enrich.ratio, avg_diff, gene.score, p_val_adj) %>%
  arrange(desc(gene.score))

  reactable(markers.tmp, defaultPageSize = 100,
          showSortable = TRUE, resizable = TRUE, highlight = TRUE, filterable = TRUE, minRows = 10,
          style = list(fontFamily = "Work Sans, sans-serif", fontSize = "12px")
  )

# saveWidget(markers.tmp, file = )

# Same genes in both comparisons (sanity check)
xtab_set <- function(A,B){
              both    <-  union(A,B)
              inA     <-  both %in% A
              inB     <-  both %in% B
              return(table(inA,inB))
            }
xtab_set(markers$nn_above_med$geme, markers$nn_below_med$gene)

Make heatmap: